Blind dispersion compensation for optical coherence tomography 
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We propose a numerical method for compensating dispersion effects in optical coherence tomography that 
does not require a priori knowledge of dispersive properties of the sample. The method is based on the gen- 
eralized autoconvolution function, and its principle of operation can be intuitively visualized using the Wigner 
distribution function formalism. 
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Since its inception about two decades ago, QJ optical co- 
herence tomography (OCT) has evolved into an increasingly 
effective and promising imaging tool, particularly in biolog- 
ical and biomedical applications. 0] The depth resolution of 
OCT is proportional to the spectral bandwidth of the light 
source as long as optical dispersion of the medium can be 
neglected. However, for increasingly broadband sources em- 
ployed to achieve ultrahigh axial resolution|3] dispersion be- 
comes relevant, broadening the coherence length of the mea- 
sured signal and thus leading to the loss of resolution. 

Generally, two strategies have been developed to combat 
dispersion effects in OCT. The first one consists in matching 
experimentally the amount of dispersion in both the arms of 
the interferometer either with a simple introduction of an ap- 
propriate compensation plate in the reference arm, 14] or by a 
skilful design of the scanning arrangement for the reference 
beam. 1 5, 6] The second strategy is to postprocess numeri- 
cally complete fringe-resolved interferograms to compensate 
for the dispersion-induced loss of resolution. |7j|8|,|9j] The nu- 
merical approach requires usually a priori knowledge of the 
dispersive properties of the medium which enter as the param- 
eters of the compensation algorithm. 

In this paper we propose a numerical method for recon- 
structing dispersion-compensated depth profiles that does not 
involve knowledge of the exact dispersive characteristics of 
the medium. This method utilizes the phase information con- 
tained in complete complex interferograms to remove delete- 
rious effects of dispersion. As we discuss later, it is equiv- 
alent to the recently demonstrated quantum OCT fTol [Till 
which uses dispersion-cancellation effect in two-photon in- 
terference. Our method offers all the advantages of quantum 
OCT, but avoids the need for non-classical light sources. The 
dispersion compensation is performed by postprocessing data 
collected in a standard OCT setup with a low-coherence light 
using a simple and straightforward in implementation numer- 
ical algorithm. 

Before passing on to the detailed discussion of the method, 
we will present its principle of operation using a simple exam- 
ple. The input for our method is the complex envelope T(t) of 
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the full analytical mutual coherence function,|7, 8, 9] param- 
eterized with the delay of the reference beam. The dispersion- 
compensated depth profile is obtained by calculating the gen- 
eralized autoconvolution of T(r), defined as: 

E w (t) = [ dT'e- 2 ™ V V(T + t')T(t - t'). (1) 



The parameter w, not related directly to the dispersion charac- 
teristics of the medium, is used to tune the performance of the 
method, and it will be typically much smaller than the light 
bandwidth. 

The autoconvolution function defined above can reveal the 
axial structure of a sample with a better resolution than the 
interferogram itself. We illustrate this with Fig.^a,b), where 
we present a reconstruction of an exemplary depth profile for 
a pair of reflective surfaces embedded in a dispersive material 
and illuminated with a broadband light. First, let us consider 
the scenario when no attempt to compensate for dispersion ef- 
fects has been made. Fig. QJa) shows the absolute value of 
the envelope |r(r)| which is the most basic way to retrieve 
the depth profile. This profile is severely broadened by dis- 
persion. For comparison we also plot in this graph a profile 
that would be obtained for a narrower bandwidth of the probe 
light, giving the optimal resolution of \T(t)\ for this specific 
medium. These two plots should be contrasted with Fig.^b) 
which shows the autoconvolution function calculated accord- 
ing to Eq. Q for several values of the parameter w. It is seen 
that the autoconvolution function reveals the location of the 
reflective surfaces as two sharp peaks with a resolution com- 
parable to the coherence length of the probe light itself. Addi- 
tionally, the convolution function contains a spurious artefact 
located half-way between the peaks. Its presence is a result of 
numerical interference of signals reflected by the two surfaces, 
but its magnitude can be quickly suppressed by increasing the 
value of the parameter w. The suppression is accompanied by 
a slight broadening of the genuine peaks in the depth profile, 
which is however significantly less severe than that affecting 

|r(r)|. 

Let us now discuss in detail properties of the autoconvolu- 
tion function R w (r) as a tool for reconstructing depth profiles. 
We consider a standard OCT setup in which broadband light, 
characterized by the power spectrum S(ui) centered around 
a frequency ujq, is split into two beams. One, signal beam 
is reflected off the sample, thus acquiring in the spectral do- 
main the response function of the sample q(oj), whereas the 
second reference beam undergoes a controlled temporal delay 
2t. Interference between these two beams yields the mutual 
coherence function, given by: 

r( T )e~ 2tuJ0T = J dujS(Lu)g(Lu)e- 2luJT . (2) 

In this formula, we separated out the phase factor rapidly os- 
cillating with the optical frequency ujq, and denoted the result- 
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FIG. 1: Reconstruction of a depth-profile of a pair of equally re- 
flective surfaces preceded by 2 cm of dispersive aqueous region 
with group velocity dispersion 15 fs 2 /mm: (a) standard interfero- 
grams |T(t)| for a coherence time T = 1.5 fim (solid line) and 
T — 5.2 ^tm (dashed line), the latter one giving optimal resul- 
tion with no dispersion compensation; (b) autoconvolution function 
H„(x) for w = 0.015 /im _1 (dotted line), w = 0.06 ^im _1 (dashed 
line), and w = 0.12 [mT 1 (solid line); (c) the chronocyclic Wigner 
function W(t, SI) of the complex interferogram envelope F(r). For 
convenience, time and frequency have been expressed in length units 
using the vacuum speed of light. The profiles in graphs (a) and (b) 
have been renormalized to the same height in order to make the com- 
parison of their widths easier. 



ing slowly-varying complex envelope as T(r). 

We will model the spectral response function g(ui) as com- 
posed of discrete contributions coming from reflective sur- 
faces within the sample characterized by reflection coeffi- 
cients r n : g(u>) — f n e 2 " p °' lJ '. The phase acquired 
by the signal field reflected from the nth surface can be ex- 
panded around the central frequency of the probe light up to 
the quadratic term: 



(3) 



ip n (uj) « Lp n (u ) + (uj - u )r n + ^(uj - uj ) 2 D n ^ 

We will incorporate the constant phase ip n (u>o) into the reflec- 
tion coefficient r n . The parameter r„ multiplying the linear 
term characterizes the position of the nth reflective surface, 
whereas D„ describes dispersion affecting the component re- 
flected from that surface. We will also assume the Gaussian 
spectrum of the probe light S(u>) oc exp(— T 2 (lo — lo^) 2 ) with 
T characterizing its coherence time. Throughout this paper, 



we will use half-width at 1/ e-maximum as a measure of the 
resolution. Within the introduced model, the complex enve- 
lope of the mutual coherence function is a sum of contribu- 
tions from the surfaces T(t) = T n (r) given by: 



r„(r) oc r„ exp 



T 2 - iD n 



(4) 



In standard OCT, the depth profile is retrieved directly from 
the interferogram as the absolute value |r n (r)|, and con- 
sequently the nth surface is visualized as a peak with a 
dispersion-broadened width \/T 2 + D 2 /T 2 . 

Let us now turn to the analysis of the information on the 
depth profile contained in the generalized autoconvolution 
function S w (t). With the complex envelope T(t) given as 
a sum of terms calculated in Eq. (0), the generalized auto- 
convolution function can be decomposed into a double sum 

^w(t) — J2mn ( r ) °f contributions obtained by insert- 

ing a product T* n (r + r')r„(r - r') into Eq. We will 
analyze separately the diagonal terms with m = n, which as 
we will see reveal positions of the reflective surfaces, and then 
the cross-terms E^ n \T) with m ^ n that are responsible 
for the artefacts in the reconstructed profile, like the one we 
have seen in Fig. \Jlb). The complete analytical expressions 
are rather complicated and we will approximate them by per- 
forming an expansion up to the leading order of w. This will 
give us an insight into relative scales of parameters involved 
in the procedure. 

For m = n, an explicit calculation yields the following 
expression for (r) in the limit when w 2 <C 1/(T 2 + 
Dl/T 2 ): 



St"")(r)(x|r„| 2 exp 



2(t - r n f 
T 2 + (D nW y 



(5) 



This expression describes a Gaussian peak located at the posi- 
tion r„ of the nth reflective surface. The width of this peak is 
given by \J\T 2 + (D„w) 2 ]/2, and in the limit when w — > 
it approaches the dispersion-free limit defined solely by the 
coherence time of the light source, equal to T/y/2. Compared 
to the standard interferogram envelope |r(r)| in the absence 
of dispersion, the peak in the autoconvolution function is nar- 
rower by a factor \/2; this narrowing is easily understandable 
as S tu (r) is quadratic in T(r). 

The purpose of introducing the parameter w is to suppress 
the cross-terms 5^7™^ (t) with m ^ n. In order to keep the 
interpretation of the mathematical expressions simple, we will 
restrict our attention to the regime when dispersion affecting 
contributions from two reflecting surfaces is comparable, i.e. 
D m w D n w D. This is the case when the bulk of dispersion 
comes from the medium preceding both the surfaces. In this 
regime, it is possible to give a simple formula for the magni- 
tude of the cross-terms: 



|H^)(r)| ex k*„rgexp - 



x exp 



2(r-f) 2 
T 2 + (Dw) 



2 ( Tm r ") 



(6) 
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where the proportionality factor is the same as in Eq. (0 and 
t = (T m + r n )/2. This formula describes a structure located 
half-way between the positions of the contributing surfaces. 
The magnitude of the structure is a function of w through the 
multiplicative factor cxp[— w 2 (r m — r„) 2 /2]. Its exponential 
dependence on w allows us to suppress efficiently the spuri- 
ous cross-terms in the autoconvolution function by setting a 
non-zero value of w, with only a slight worsening of the reso- 
lution exhibited in Eq. (0. The supression of the cross-terms 
requires \w\ exceeding l/|r m — r n \ and therefore is more ef- 
ficient for larger separation between the peaks. 

The operation of the blind dispersion compensation method 
can be understood intuitively with the help of the chronocyclic 
Wigner distribution function,! 12] defined for the complex in- 
terferogram envelope T(t) in the standard way as: 

W{t, ®) = \ I dT ' e- 2inr 'T*{T + t')T{t - t') (7) 

In Fig. [flc) we depict the Wigner function for the example 
discussed in this paper. The Wigner function contains two 
peaks corresponding to the two reflective surfaces, and an 
oscillating interference pattern located half-way between the 
peaks. This pattern is a signature of coherence between the 
two reflections. Dispersion introduces time-frequency cor- 
relations which result in a tilt clearly seen in Fig. ^ c )- In 
standard OCT the depth profile is retrieved as |r(r)| which is 
given as a square root of W(t, 57) integrated along the ft axis. 
Then the dispersion-induced tilt severely deteriorates the res- 
olution. However, taking a cross-section through the Wigner 
function along a horizontal line for O = 0, yields a dispersion- 
free profile. The problem with this profile is that it contains 



strong contributions from the interference pattern! 13] which 
were washed out in |r(r)| due to the integration over the fre- 
quency variable. The answer to this problem is to carefully 
average the profile over a range of frequencies ft. This is ex- 
actly the purpose of the generalized autoconvolution function 
in Eq. Q, which can be rewritten in the Wigner formalism as: 

Z w (t) = J dfi W(t, Q)e^ 2 / 2w2 . (8) 

The averaging along the frequency axis over an interval de- 
fined by the parameter w rapidly washes out the contribution 
from the interference pattern, while nearly retaining the width 
of the genuine peaks. The width of the interval must be larger 
than the spacing of the interference pattern, which in turn is 
inversely proportional to the separation between the peaks. As 
shown by Abouraddy et al. | 10], the effects of dispersion can 
be removed using a quantum effect of two-photon interfer- 
ence. Interestingly, the depth profile obtained from the joint 
detection of two photons in their case yields the same order 
correlation function as E w (t). This equivalence is related to 
the fact that in the quantum OCT scheme only one of the two 
photons passes through the dispersive medium. 

In conclusion, we have shown that the generalized autocon- 
volution function calculated from the complex interferogram 
envelope can reveal location of reflective surfaces in a disper- 
sive medium with a resolution reaching the coherence length 
of the employed light itself. We expect that in the case of 
more complex depth profiles, the oscillatory character of the 
artefacts will lead to similar suppression efficiency as that in 
the simple numerical example discussed here. 
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